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Current reduced-order thermal model for cryogenic propellant tanks is based on 
correlations built for flat plates collected in the 1950’s. The use of these correlations suffers 
from inaccurate geometry representation; inaccurate gravity orientation; ambiguous length 
scale; and lack of detailed validation. This study uses first-principles based CFD methodology 
to compute heat transfer from the tank wall to the cryogenic fluids and extracts and correlates 
the equivalent heat transfer coefficient to support reduced -order thermal model. The CFD tool 
was first validated against available experimental data and commonly used correlations for 
natural convection along a vertically heated wall. Good agreements between the present 
prediction and experimental data have been found for flows in laminar as well turbulent 
regimes. The convective heat transfer between the tank wall and cryogenic propellant, and 
that between the tank wall and ullage gas were then simulated. The results showed that the 
commonly used heat transfer correlations for either vertical or horizontal plate over-predict 
heat transfer rate for the cryogenic tank, in some cases by as much as one order of magnitude. 
A characteristic length scale has been defined that can correlate all heat transfer coefficients 
for different fill levels into a single curve. This curve can be used for the reduced-order heat 
transfer model analysis. 


Nomenclature 


M I. Introduction 

odeling propellant tank dynamics including slosh, pressurization, propellant conditions, and 
thermodynamics of heat and mass transfer across gas-liquid interface for liquid propellant based space 
vehicles is very important. In the development of Saturn stages, such as S-IC, S-II, S-IVB, in the 1960s, 
the analyses were done typically using reduced thermal models [1,2]. With the advancement of Computational Fluid 
Dynamics (CFD) technology, it is possible to model the propellant tank dynamics from a first-principles based 
approach. The reduced thermal models [1-4] have the advantage of fast turnaround time for parametric studies and 
are easy to use. Recently, an improved analytical thermal modeling for propellant tank has been developed at NASA 
MSFC [4]. The model validation has been proven successful against previous models [3] and various Saturn era test 
data and reasonably successful against more recent LH2 tank self-pressurization ground test data. 
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The reduced thermal models for cryogenic tanks [2,4] typically divide the propellant tank into nodes and 
implement an empirical correlation for heat transfer between different nodes. For example in the model of Reference 
[4], fives nodes are used as shown in Figure 1, with a tank wall node exposed to ullage gas, an ullage gas node, a 
saturated propellant vapor node at the liquid-vapor interface, a liquid node, and a tank wall node exposed to liquid. 
The conservation equations for mass and energy are then applied across all the node boundaries and, with the use of 
perfect gas assumptions, explicit solutions for ullage and liquid conditions are derived. 
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Figure 1: Diagram of control volume for the analytical thermal model [4]. 


From a fluid physics point of view, the heat transfer between the tank wall node exposed to the ullage gas 
and the ullage gas node, and that between the liquid node and the tank wall node exposed to the liquid, are all by 
natural convection which involves complex fluid dynamics. In the reduced thermal model, the calculation of heat 
and mass transfer between the cryogenic propellant tank wall and the internal fluids was based on correlation 
dependencies created from experimental data. Some of these correlations were built long ago [2] (condensation on 
the vertical wall - 1960) and some more recently; but some of them still need to be constructed. One of the most 
important parameters is the correlation for heat transfer from the tank dome to the ullage gas and from the tank wall 
to the liquid (see Fig. 1). At the present time, calculations are made by converting spherical or elliptical domes to 
equivalent flat plates and using correlations for the flat plates. Even so, for the upper domes or for the downward- 
facing plates, the correlations don’t account for the presence of a cylindrical barrel section, which must be 
considered in practical design. New correlations are still required to properly account for the heat transfer between 
the tank wall and the propellants inside the tank from which a reliably reduced model can be formulated. 


The application of existing correlations developed from flat plates to cryogenic tanks suffers from: 

• Inaccurate geometry representation (spherical and cylindrical surfaces are approximated by a flat surface); 

• Inaccurate gravity orientation (an inclined surface is approximated as either a vertical or a horizontal 
surface); 

• Ambiguous length scale (tank diameter, barrel height, or tank height); 

• Ambiguous computation domain as the correlated heat transfer coefficients are constructed for a flat 
surface in an infinite medium, whereas the application is for the finite space inside a tank; and 


• Lack of detailed validation (the temperature field inside fluid is stratified rather than uniform as assumed in 
the reduced-order model, and validation must be made with respect to volume averaged temperature which 
is difficult to obtain from experiment) 

The objective of this study is to use a first-principles based CFD technique to compute heat transfer from the 
tank wall to the cryogenic fluids, and extract and correlate the equivalent heat transfer coefficient to support 
reduced-order thermal model development. There are many advantages of the current approach: 

• ability to use the exact tank geometry; 

• a validated CFD model; 

• physics based length scale; and 

• a heat transfer data extraction procedure that is consistent with the reduced-order thermal model to ensure 
energy conservation. 


II. Computational Fluid Dynamics (CFD) Solution 

The basic CFD software used to study the heat transfer inside the cryogenic tank is the commercial code, CFD- 
ACE+, which was originally developed by CFD Research Corp., and is currently owned and distributed by ESI 
Group USA. CFD-ACE+ solves conservation of mass, momentum (Navier- Stokes equations) and energy. The mass 
and momentum continuity equations can be generally written as: 
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where p is the density, u, v, and w are the velocity components in x, y, and z directions. V is the velocity vector, p 
the static pressure, t the time, and Smx, Sm y and Sm z are the momentum source in x, y, and z directions, respectively. 


For Newtonian flows, the viscous stresses (lij) are proportional to the deformation rates of the fluid 
element. The nine viscous stress components can be related to velocity gradients to produce the following shear 
stress terms: 
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For the natural convection problem under this study, the momentum source includes the buoyancy force 
applied through the Boussinesq approximation. 

For heat transfer problems, CFD-ACE+ solves the total enthalpy equation derived from energy 
conservation: 
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where h 0 is the total enthalpy and is defined as: 
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where i is the internal energy and is a function of the state variables of density p and temperature T. K e ff is the 
effective thermal conductivity of the material. In laminar flow, K e ff is the thermal conductivity of the fluid, K. In 
turbulent flows: 


K 
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with ju t being the turbulent viscosity and <j t the turbulent Prandtl number. 
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The most relevant fluid dynamics features of CFD-ACE+ are: 

• Pressure -based, strongly implicit algorithm solving Navier-Stokes equations for all flow speeds; 

• High order TVD scheme and high order time-stepping for high spatial and temporal resolution; 

• Unstructured, mixed element grids for complex geometrical problems, including tetra-, hexa-, adaptive 
Cartesians, pyramids, prisms, polyhedral and cells with multiple hanging nodes; 

• Co-located (non-staggered) grid topology; 

• A space-conserving moving grid formulation that allows the treatment of moving -deforming grid systems; 

• Comprehensive set of boundary conditions; 

• Stationary, deforming, moving, rotating frame of reference; 

• Various turbulence models for Reynolds Averaged Navier-Stokes Equation, including: 

■ Standard k-s Model of Launder and Spalding [5] 

■ RNG k-s Model 

■ Realizable k-s (RKE) Model 

■ Kato-Launder k-s Model 

■ Low Reynolds Number k-s Model (Chien) 

■ Two-Layer k-s Model of Chen & Patel [6] 

■ k-CD Model 

■ k-coSST Model 

■ Spalart-Allmaras Model 

• Various subgrid-scale stress models for Large Eddy Simulations 

■ Smagorinsky SGS model 

■ Germano’s Dynamic Sub grid- Scaled model 

■ Menon’s Localized Dynamic Subgrid-Scale Model 


III. Validation Study of Natural Convection over a Vertical Plate 

Before we apply CFD technology to the cryogenic tank heat transfer problem, it is very important to make a 
systematic validation study. As the empirical correlation of heat transfer coefficient for the flat plate has been widely 
used in the reduced thermal model, it is crucial to verify that CFD can reproduce the classical experimental data. The 
first task here is to perform validation for natural convection along a vertical plate in both laminar and turbulent 
regimes with air as the working medium. 



A. Validation Model and Fluid Physics 

Figure 2 shows the fluid dynamics process of natural convection along a vertical plate. Here a fluid medium is 
subjected to a warm wall temperature. The warm wall heats the near- wall fluid. Near- wall fluid flows in +x direction 
due to buoyancy effects. As the fluid flows along the wall, it draws in more heat (T is a function of y). Thermal 
boundary layer forms and grows thicker because of local heat transfer q”(x). More fluid is drawn in from the fluid 
medium, and the process continues. 


There are several important non-dimensional parameters for the natural convection problem. 

a. Prandtl Number. Prandtl number is the ratio of viscous diffusivity, v, to thermal diffusivity, a, and is 
defined as: 
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Rayleigh Number. Rayleigh number is the ratio of buoyant forces to viscous forces multiplied by the ratio 
of viscous diffusivity to thermal diffusivity. Rayleigh number is defined as: 
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where L is the characteristic length, p is the thermal expansion coefficient, and AT is the temperature 
difference. 


c. Local Nusselt number. Local Nusselt number is the ratio of convective to conductive heat transfer. The 
typical formulation for Nusselt number is Nu = f(Pr, Ra). Nusselt number is given by: 


Nu = 


hx 
K ’ 
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Figure 2: Fluid physics of natural convection over a vertical flat plate. 


B. Previous Experimental Data and Empirical Correlations 

For natural convection over a vertical flat plate, several experiment-derived datasets exist correlating the non- 
dimensional Nusselt number as a function of Ra as shown in Figure 3. During the early development of an analytical 
similarity solution for natural convection along a vertical plate, these data were used to validate the asymptotic 
analytical solution [7]. For example, the boundary layer similarity solution gives 1/4 power law for laminar flow, 
and 1/3 power law for turbulent flow, as seen in Figure 3. These empirical correlations have been widely used in the 
reduced thermal model analysis for cryogenic tank thermal analysis [2, 4]. The general form of the averaged natural 
convection expression is: 


Nu = C(Ra) n . 


( 12 ) 


Based on the orientation of the flat plate to the gravity vector, and Ra number, C, and n take on different 


values. 



Figure 3: Experimental data and empirical correlations for natural convection over a vertical flat plate from 

surveyed literature [7]. 


Experimental Data Correlations 

For the vertical surfaces: 

For Laminar Range: 10 4 <Ra<10 9 ; C=0.59, n=l/4. 

For turbulent Range: 10 9 <Ra<10 13 ; C=0.13, n=l/3. 

For the horizontal surfaces: above hotter surface, below cooler surfaces 
For Laminar Range: 10 5 <Ra < 2xl0 7 ; C=0.54, n=l/4 

For turbulent Range: 2xl0 7 < Ra < 3xl0 10 ; C=0.14, n=l/3 

For the horizontal surfaces: below hotter surface, above cooler surfaces 
For Laminar Range: 3xl0 5 <Ra < 3xl0 10 ; C=0.27, n=l/4 
No correlation in the turbulent regime. 

The above experimental data correlations will be used for the present validation study. 


C. CFD Model and Boundary Conditions 

The CFD simulation domain is shown in Figure 4. As the experimental data were collected in an infinite 
medium, the size of the domain is selected such that the far field boundary condition has limited impact on the 
results. A 2D model was built with an assumption that the plate is wide enough such that the edge effect can be 
neglected. The flat plate itself is assumed to be infinitely thin so that only half of the domain needs to be considered. 
The condition on the flat plate is a constant temperature and non-slip velocity. Below the plate, it has a vertical 
symmetry condition. At far field down below the plate, the pressure is constant, and the temperature is the 
freestream value. On the right-hand side of the domain, a symmetry condition is imposed. 



Fixed Pressure 

Figure 4: CFD simulation model for heat transfer over a vertical flat plate. 

The computational grid is shown in Figure 5, where the mesh is packed near the plate. There is a total of 90k 
cells in the baseline grid. 
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Figure 5: Computational grid for the vertical flat plate model. 


D. Validation for Laminar Regime 

For the validation studies of laminar flow over a flat plate, we consider Ra number up to 10 6 . From the equation 
(10), one can see that Rayleigh number is a function of gravity acceleration g. Under the same thermal conditions, 
different gravity levels can give different Rayleigh numbers. By changing gravity input, one can change the 
Rayleigh number. For air, the properties are given in Table 1 under ambient temperature of T=300°K and ambient 
pressure of P= 1.0x1 0 5 N/m 2 . 


Table 1. Air Properties Used in the Validation Study. 


Property 

Value 

Density, p (kg/m3) 

1.225 

Viscosity, p (kg/m*s) 

1.79E-5 

Specific Heat, Cp (J/kg*K) 

1005 

Thermal Conductivity, K (W/m*K) 

.0257 

Thermal Expansion Coefficient, (3 (1/K) 

.00343 




( 13 ) 


To extract the area averaged heat transfer rate, the local heat flux into the wall is first calculated as: 


q'\x') = K f 


dT 

~dy 


The mean heat flux up to the point x is the result of averaging the above heat flux: 


q'\x) = 



Based on the mean heat flux, one can compute the mean Nusselt number as: 


Nu(x) = 


q" x 
( X~TJK f 
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The length scale in the Rayleigh number is based on the above x. 

Given in Figure 6 is the local heat transfer computed using equation (14) after the CFD solution has reached 
steady state. Both 1 st and 2 nd spatial schemes are applied. Near the leading edge of the plate, the boundary layer is 
very thin which gives rise to a locally high heat flux. After the leading edge, the heat flux is nearly uniform. 

By integrating the surface heat flux value as the above equation (14) and using the definition of (15), the present 
CFD simulation of heat transfer rate in terms of Nusselt number as a function of Rayleigh number is given in Figure 
7. The experimental data and empirical correlation based on thermal boundary layer solution [7] is also shown in the 
same figure for this laminar flow regime. One can see that from Ra=10 4 and up, all three agree with each other very 
well in term of heat transfer value. Also, all three show the same trend, % power of Rayleigh number. This implies 
the laminar boundary layer as the dominant heat transfer mechanism. When Rayleigh number is below 10 4 , one 
starts to see the deviation of empirical correlation relative to the experimental data and the present CFD simulation. 
It is expected that at low Rayleigh number, the heat transfer mechanism is conduction dominated, and Nusselt 
number should approach a constant. That is indeed what experimental measurements and the CFD predictions show. 
In terms of the present CFD simulation, one can see fair good agreement with experimental data even down to 
Rayleigh number order of 1 . 
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Figure 6: Local heat flux from the present simulation using different differencing schemes. 



Figure 7: Comparison of mean Nusselt number from the present CFD prediction, experimental data, and 

empirical correlations. 

To ensure the solution is mesh- independent, the simulations are made with additional two finer meshes of 120k 
and 240k. The results are displayed in Figure 8. It is clear that the present baseline 90k cell gives a grid-independent 
solution. 



Figure 8: Mesh independent study for laminar natural convection over a vertical flat plate. 

E. Validation for Turbulent Regime 

When Rayleigh is above 10 9 , flow is turbulent. To model the enhanced heat transfer due to turbulence, one can 
use one of many turbulence models in CFD-ACE+. The following four turbulence models have been assessed 
against experimental data: 

• Standard k-s model of Launder and Spalding [5]; 

• K-cd model; 


• STT model, which blends k-s near wall and k-co away from the wall; 

• two layer k-s model of Chen & Patel [6], which blends high Reynolds number and low Reynolds number 
model. 


The results are shown in Figure 9. It can be concluded that: 

• k-s model diverges from all other models and experimental data at high Ra; 

• SST model is highly over-predictive at low Rayleigh number; 

• two layer k-s and k-co models are very competitive. 


Based on the above observations, it was decided that the two-layer k-s model will be used for the subsequent 
simulations, as the model is designed for resolving viscous and thermal boundary layers without using wall 
functions. Another feature of two layer k-s is that it can resolve the effect of surface curvature and flow separation, 
which could occur for cryogenic tank problems. 



Figure 9: Evaluation of turbulence models for natural convection over a vertical flat plate. 


Figure 10 shows the comparison of heat transfer in the turbulent flow regime using the two layer k-s model, 
along with the experimental data and empirical correlations. One observes that the CFD simulations match very well 
with the experimental data from low Rayleigh number of Ra=10 4 up to high Rayleigh number of Ra=10 n . The 
empirical correlation of heat transfer with a 1/3 power to Ra in high Rayleigh number regime and % power to Ra in 
lower Rayleigh number regime are all displayed from the simulation results. There is some scatter in the 
experimental data in the transition regime, but all the essential physics of heat transfer are well represented in the 
present numerical study. This builds a sound foundation for the heat transfer study inside cryogenic tanks. 



Figure 10: Validation of two layer k-e model for natural convection over a vertical plate. 

Figure 1 1 shows the mesh-independence study for turbulence simulation. Indeed solutions for different meshes 
show little or no difference, and mesh independence is achieved at the 90,000 cell level. 



Figure 11: Heat transfer sensitivity to grid mesh sizes for the two-layer k-e model. 

IV. CFD Extraction of Heat Transfer Coefficient 

Analysis of heat transfer inside cryogenic tanks is crucial for the design of space launch vehicles. There is, 
however, a lack of such experimental heat transfer data available for cryogenic tanks. Our current approach is to use 
a validated CFD tool to predict and to correlate the natural convection heat transfer between the liquid phase and 



the tank wall, and between gas phase and the tank wall. The successful validation exercise outlined in the 
previous sections provides a sound foundation for such predictions. 


A. Technical Approach for Heat Transfer in Cryogenic Tank 

The general approach is to partition the tank into several sub -models, and use CFD to extract the heat transfer 
coefficient for each sub-model. The extracted heat transfer coefficient will then be employed in a reduced thermal 
model analysis. 

Before extracting the heat transfer coefficient, it is important to review the reduced thermal model used in the 
cryogenic tank analysis [2,4]. The typical reduced thermal model divides the tank into 5 nodes as shown in Figure 
11 . 

Tank wall node exposed to ullage gas; 

Ullage gas node; 

Saturated propellant vapor node at liquid-vapor interface; 

Liquid node; and 
Tank wall node exposed to liquid. 
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Figure 12: Reduced thermal model for cryogenic tank. 

The energy conservation equation for the ullage gas node (assuming no mass transfer) can be written as: 


d (Pc c P G T a) _ 

dt ^ wg ~ gs ' 


= Q -Q 


(16) 


where Q wg is the heat flux from the tank wall to the ullage, and Q gs is the heat flux from the ullage to gas-liquid 
interface. And Q wg and Q gs are given by: 


Qwg h wg A wg AT , 
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The energy conservation for the liquid node is 

S(Pl c p l T l) 


dt 
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where Q w i is the heat flux from the tank wall to liquid, and Qi s is the heat flux from liquid to gas -liquid interface. 
And Qwi and Ql s are given by: 


Qwi — h v/l A wl !S3 , 

(20) 

Qis = MAI . 

(21) 


Our goal is to determine the heat transfer coefficients of h wg , h gs , h w i, and hi s . 

B. Extraction of Heat Transfer Rate by Steady State Simulation 

Our experience indicated that the extraction of a converged heat transfer coefficient from a transient simulation 
is difficult, and a meaningful extraction of heat transfer rate between tank walls and the internal fluids is by a steady- 
state approach. A steady-state solution to heat transfer problem requires a heat source and a heat sink. In the current 
cryogenic tank problem, the tank wall can be taken as a heat source. The question is where the heat sink is located. It 
is noted that as the liquid phase has a much higher heat capacity, and the gas-liquid interface is maintained at the 
saturated temperature, it is reasonable to assume that gas-liquid interface behaves as the heat sink as illustrated in 
Figure 13. With high and low temperatures specified, the ullage gas is the media between the tank wall and gas- 
liquid interface to convect the heat. 


Isothermal tank wall "^WG 



^tank 

Figure 13: Model for steady solution of tank heat transfer problem. 

The total heat transfer between the tank wall and the ullage gas can be expressed as: 

Q wg = KgKgiTwc ~ t gl) ■ 

The total heat transfer from ullage to the gas-liquid interface is: 


( 22 ) 


Q„=KA t ,(T wa -T aL ). (23) 

At the steady state, the above two fluxes are equal, so one has the relationship of: 

U A — h A 

n wg jr Kvg n gs rx gs' ( 24 ) 

Once heat transfer coefficient from the tank wall to ullage gas h wg is determined, the heat transfer from ullage gas to 
the gas-liquid interface can also be determined be as: 



(25) 
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A similar procedure can be applied to heat transfer coefficient between the tank wall and liquid phase h w i 
and that between the liquid phase and gas-liquid interface hi s . 

The extraction procedure is as follows: 


1. 

2 . 


3. 


4. 


5. 


6 . 


Solve steady-state heat transfer problem with constant tank wall temperature (Twg) and constant gas-liquid 
surface temperature (Tgl); 

Compute local heat flux to the wall: 

4 wall = ^eff (26) 

Compute the mean heat transfer coefficient based on the total exposed tank surface area: 


h = 


J q wa n dA 
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(27) 
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Define an equivalent total heat conduction flux based on a characteristic length L and fluid conductivity Kf as: 


Qcond ~ K 


(t wg -t gl ) a 


L W8 

Compute the non-dimensionalized Nusselt Number as 

Qconvect hL 

Qcond ^ f 

Correlate Nusselt number as a function of Rayleigh number 


Nu = 


(28) 

(29) 


Using the above approach, the natural convection in the gas phase and in the liquid phase can be separated. 
A generic tank with equal sized upper and lower domes is considered as illustrated in Figure 14. Here the barrel 
(cylindrical) section has a height of tank radius R. Three liquid level fill positions are considered: a) fill level to top 
of the barrel section, where the gas phase fills the upper dome (Figure 14a); b) fill level to the half-height of the 
barrel section, so that the liquid and gas phases have the equal volumes (Figure 14b); and c) fill level to the lower 
dome only, which is opposite to case a) (Figure 14c). The problem can now be separated into solving natural 
convection in the liquid, geometry corresponding to red color and natural convection in the gas phase with a gray 
colored geometry in Figure 14. 








(a) (b) (c) 

Figure 14: Three fill levels and models under consideration for natural convection inside a propellant tank. 


V. Natural Convection of Liquid Propellant inside Tank 
A. Simulation Model 

The computational models for natural convection inside the liquid phase at the three fill levels shown in Figure 
14 are given in Figure 15. The flow and temperature are assumed to be axis-symmetric so that a 2D model is 
sufficient. The grid is packed near wall with a total cell number varying from 28K to 73 K. A symmetry boundary 
condition is applied along the center axis. The gas-liquid interface is specified at a constant temperature of 21.46°K, 
while the tank wall is at a constant temperature of 30°K. The physical properties are assumed to be constant and are 
taken as the values at the reference temperature of 21.46°K. The buoyancy force is applied in the momentum 
equation through the Boussinesq approximation, which implies that the density differences are sufficiently small to 
be neglected, except where they appear in the terms multiplied by g, the acceleration due to gravity. 



Figure 15: CFD model for the liquid convection inside tank with different fill levels. 




The no-slip boundary condition is enforced on the gas-liquid interface as well as on the wall. In the reality, 
the gas-liquid interface satisfies a continuous shear stress condition along the interface which falls in between no- 
slip and slip conditions. From that perspective, it will require the solution of both gas and liquid phases together. 
Since under steady state condition, the heat flux from tank wall to the liquid is the same as heat flux from liquid to 
the gas-liquid interface, the use of a no-slip boundary condition at interface will only affect the local temperature 
distribution. It has a small effect on the total heat flux, as we are more interested in the heat flux near the tank wall. 
For simplification, a no-slip condition at the interface is justified. 

Table 2 lists the physical properties used in the simulation. As evident from Equation (10), many 
parameters can affect the value of Rayleigh number. For example, the characteristic length, L, which is typically the 
tank radius, the temperature differential (AT), fluid properties, and gravity condition (g). Under different durations 
of a flight, the total acceleration can be very different, varying from 1 g on the ground to 4g during the ascent of the 
flight, and to micro-g in space. In order to cover a wide of Rayleigh number, the gravity value is varied in the 
simulation. 

T able 2. Liquid Properties under Consideration for the Natural Convection. 


Property 

Value 

Density, p (kg/m3) 

69.565 

Viscosity, p (kg/m*s) 

1.2228xl0" 5 

Specific Heat, Cp (KJ/kg*K) 

10.419 

Thermal Conductivity, K (W/m*K) 

0.10491 

Thermal Expansion Coefficient, (3 (1/K) 

0.01754 

Reference Temperature (K) 

21.46 

Gravity g (m/s 2 ) 

2.79xl0" 9 to 9.8xl0 2 


B. Natural Convection for Liquid Filled to Lower Dome 

First we consider that case when the lower dome is filled with liquid, Case c of Figure 14. As Rayleigh number 
is a non-dimensional parameter, the Rayleigh number is allowed to vary from 10 7 to 3.0x1 0 9 after scaling using the 
radius of these tanks. The temperature fields and streamlines for these tanks are given in Figure 16. At Rayleigh 
number of 1 0 7 , it is observed that a thermal boundary develops along the tank wall starting from the center line all 
the way to the gas-liquid interface. The fluid temperature increases by absorbing heat from the wall and the fluid 
flows up due to buoyancy to the gas-liquid interface. After hitting the cooler gas-liquid interface, the fluid gives up 
heat and is cooled while returning back along the center of the tank. Due to axisymmetry, the streamlines of natural 
convection are in the form of a torus. The thermal boundary layer becomes thinner with the increase of Rayleigh 
number for Ra=3xl0 7 and 3x1 0 9 . Based on the above study of a flat plate, two-layer turbulence model is applied 
when the Rayleigh number is above 10 8 . Numerical experiments show an independence of the solution on the use of 
turbulence model for Ra below 10 8 . Figure 17 displays the temperature fields and streamlines under even higher 
Rayleigh number. The thermal boundary layer is so thin that it is almost invisible. It should be noted that the grid is 
very fine near the tank wall, and there are at least 5 points inside the thermal boundary layer. In general the flow 
field is about the same showing a torus flow pattern. 

The distribution of local heat flux based on Eq. (26) is plotted in Figure 18. A singularity is seen as the heat flux 
near the location where tank wall meets the gas -liquid interface. At the center line portion of the tank wall, a higher 
heat transfer is noticed due to the direct impingement of cold flow to the wall. The increase of thermal boundary 
thickness is observed from the temperature field, which leads to the reduction of heat flux until a sudden spike at the 
singularity point. 
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(c) Ra=3x10 9 , Two Layer Turbulence Model 

Figure 16: Temperature field and streamlines for natural convection inside liquid filled to the lower dome 

under laminar Rayleigh numbers. 



(a) RasIxlO 13 , Turbulence Model 
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(b) Ra=3x10 13 , Turbulence Model 
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(c) Ra=10 15 , Turbulence Model 

Figure 17: Temperature fields and streamlines for natural convection inside liquid filled to the lower dome at 

high Rayleigh numbers. 
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Figure 18: Distribution of local heat flux along the tank wall and the related temperature field. 
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With the computed heat flux, it is now possible to extract Nusselt number based on the procedure discussed 
earlier. As there is no direct length scale, the tank radius will be used in the calculation of Rayleigh and Nusselt 
numbers. The consolidated heat transfer rate is shown in Figure 19. The experimental correlation for the vertical 
plate under laminar and turbulent conditions and for the horizontal plate heated below are also shown in the same 
figure. 


Extracted Nusselt Number 
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Figure 19: Nusselt number as a function of Rayleigh number for the lower dome filled with liquid. 


It is observed that: 

1. The application of the empirical curve correlated from vertically heated plate to a tank dome, as 
practiced in reduced order analysis, over-predicts heat transfer rate by as much as 5 times. 

2. Empirical correlation for a surface heated from below fits the lower dome heat transfer characteristics 
very well. But the correlation is only available for laminar regime (Ra < 3x1 0 10 ). 

3. CFD predictions are able to provide a reliable heat transfer coefficient for liquid propellant tank. 


C. Natural Convection for Liquid Filled to Half Barrel Height 

The Nusselt number based on the tank radius is plotted in Figure 20 at selected Rayleigh numbers, which are 
also based on tank radius. It is seen that the use of empirical correlations for either vertical plate or horizontal plate 



will over-predict heat transfer rate. The heat transfer rates from the present CFD predictions for the liquid fill levels 
at lower dome and at half radius height in barrel section don’t collapse into a single curve using tank radius as the 
characteristic length. 


Extracted Nusselt Number 
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Figure 20: Nusselt number as a function of Rayleigh number for liquid level to 1/2R in barrel 

section. 

D. Natural Convection for Liquid Filled to Full Barrel Height 

The Nusselt number based on the tank radius is plotted in Figure 21 at selected Rayleigh numbers, which are 
also based on tank radius as the length scale. It is seen that the use of empirical correlations for either vertical plate 
or horizontal plate tends to over-predict the heat transfer rate. Surprisingly, even under this condition with much 
more exposed vertical surface in the tank, heat transfer rate very much follows the trend of a horizontally heated 
plate rather than that of a vertically heated plate. One also notices that the CFD-predicted heat transfer rate for all 
three fill levels: up to lower dome, at half radius height and full radius height in barrel section, do not collapse into 
single curve using radius of tank as the characteristic length. 
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Figure 21: Nusselt number as a function of Rayleigh number for liquid level filled to a radius height in barrel 

section. 


E. Correlation Development for Natural Convection of Liquid Propellant inside Tank 

Based on the above study, we find that using the tank radius as the characteristic length, the resultant Nusselt 
number has different curves for different fill levels. Without a generalized correlation, it is difficult to apply the 
present results in engineering applications. The challenge is to find a characteristic length that will collapse the 
curves into a single correlation. In natural convection with a boundary layer, such as the case of flat plate, the 
vertical height is one of the important length scales. Here we propose a representative length scale for different fill 
levels as shown in Figure 22. This length scale is the wetted vertical height inside the tank. With this characteristic 
length, we have the following conversion for Nusselt and Rayleigh numbers: 


Nu l 



(30) 


f j Y 


and 


Ra L = 


\Rj 


Ra, 


(31) 


Using the above definition, the resultant Nusselt number, and Rayleigh number for different fill levels are 
plotted in Figure 23. It is clearly seen that all the data collapse into a single curve. This curve can be fitted into the 
following correlations: 



Figure 22: Representative length scale for different fill levels in the liquid tank. 
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Figure 23: Curve fit using wetted liquid height as the characteristic length. 


Nu l = 0.642 (Ra L f 6 Ra L < 10 7 

Nu l = 0.167(i?a L ) 1/4 10 7 < Ra L < 10 10 (32) 

Nu l = 0.00053 (Ra L ) l/2 10 10 < Ra L < 5xl0 13 

The above correlation is consistent with the previous study of a vertical flat plate in that in the laminar regime, 
the Nusselt number has a % power to the Rayleigh number as predicted by analytical solution for laminar flow using 
similarity solution. In the turbulence regime, the Nusselt number grows faster inside the tank than that for the case of 
flat plate in an infinite medium. The current correlation displays a power of l A rather than 1/3 with Rayleigh number. 
In the conduction regime for Ra < 10 7 , it has a power of 1/6. 

The present correlation is compared with the previous experimental correlations for flat plate in Figure 24. One 
notices that the present curve trends closely with the heat transfer correlation for the horizontal plate geometry 
heated from below. However, the new correlation developed from this CFD study significantly increases the 
accuracy of heat transfer prediction and additionally shed light on the use of the liquid height as a viable 
characteristic length for non-dimensionalization. 
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Figure 24: Present correlation for natural convection of liquid propellant inside a tank compared to that for 

heated flat plate. 


VII. Natural Convection of Ullage Gas inside Tank 

A. Simulation Model 

Next, we consider the natural convection of ullage gas above the gas-liquid interface. Here the tank wall is at a 
higher temperature and is above the gas-liquid interface as illustrated in Figure 14. By changing the gravity vector, 
the same three sets of grid for liquid phase can be used for the ullage gas phase. The computational model for the 
three fill levels is illustrated in Figure 25. The flow is assumed to be axisymmetric so that a 2D model is sufficient. 
The tank has the same size as the one used for liquid phase convection with a radius of 33.45 inches (0.84693m). 
The grid is packed near the wall with a total cell number varying from 28K to 73 K. A symmetry boundary condition 
is applied along the center axis. The gas-liquid interface is specified at a constant temperature of 120°K, while the 
tank wall is at a constant temperature of 140°K. The gas-liquid interface temperature in the present gas phase 
convection is set at a different temperature from the case of the liquid phase convection. This is to reflect the fact 
that the ullage gas is typically stratified and has a higher temperature than that of the liquid. It is should be noted that 
under the assumption of constant physical properties and Boussinesq approximation, heat transfer rate is only a 
function of non-dimensional Rayleigh number. The selection of the above temperatures is from typical application 
encountered in cryogenic fluids. The physical properties are assumed to be constant and are taken at 130°K. The 
buoyancy force is applied in the momentum equation through Boussinesq approximation. The no -slip boundary 
condition is enforced at gas -liquid interface as well as on the wall. 

Table 3 lists the physical properties of ullage gas in the present simulations. It should be noted that to cover the 
whole range of Rayleigh number, from 10 4 to 10 12 , the gravity value is different from that in the liquid case, as the 
gas physical properties are very much different from those of liquid. 



Table 3. Ullage Gas Properties under Consideration for the Natural Convection. 


Property 

Value 

Density, p (kg/m3) 

0.50895 

Viscosity, p (kg/m*s) 

1.138xl0" 5 

Specific Heat, Cp (KJ/kg*K) 

0.5193 

Thermal Conductivity, K (W/m*K) 

0.088 

Thermal Expansion Coefficient, P (1/K) 

0.01535 

Reference Temperature (K) 

130 

Gravity g (m/s 2 ) 

7.9xl0" 5 to 7.9xl0 3 





Figure 25: CFD model for the ullage gas convection inside tank with different liquid fill levels. 


B. Natural Convection of Ullage Gas 

First we consider the case when the liquid is filled to the bottom of the upper dome. The computational domain 
is simply an upper dome. At the bottom of the upper dome is the gas-liquid interface with a constant temperature of 
130°K. Along the tank wall, the temperature is fixed at a higher temperature of 140°K. Figure 26 shows the 
temperature fields and streamlines for different Rayleigh numbers. At Ra=10 4 , the temperature field is stratified 
vertically with low convection current along the wall. As the heated tank wall is mainly situated on the top of the 
cold gas-liquid interface, the lighter fluid (hotter gas) is on the top of the heavier fluid (colder gas), so the 
convection is expected to be insignificant. That is indeed the case as one observes from the streamline distribution. 
The temperature field is essentially conduction dominated with convection arising near the dome -cylindrical section 
interface where the temperature gradient becomes normal to the gravity vector (unstable condition). One sees 
isothermal contours stratified except near the corner where tank wall meets the gas -liquid interface. As Ra increases 
to 10 8 , the temperature field is again stratified but with a higher gradient implying an increase in the total heat 
transfer. At Ra reaches a value of 10 10 as shown in Figure 26(c), cells (similar to Benard cells) are visible from the 
streamlines. It is expected that heat transfer will be enhanced due to the convection of these small cells on the 
interface. 
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(c) Ra=10 10 , Turbulence Model Bernard Cells 

Figure 26: Temperature fields and streamlines for natural convection of ullage gas inside a tank with liquid 

filled to the bottom of the upper dome. 


Figure 27 shows the temperature and flow fields for the liquid level filled to half and full radius heights. At 
Ra=10 4 , the temperature field is stratified and shows a conduction pattern. However, due to the difference in the 
barrel height, the heat transfer rate for the two cases is expected to be different. 
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(b) Ra=10 4 , Laminar 

Figure 27: Temperature fields and streamlines for natural convection of ullage gas inside a tank with liquid 

filled to half of radius height in barrel section. 




C. Correlation Development for Natural Convection of Ullage Gas inside Tank 

The heat transfer coefficients are extracted from the above three liquid fill levels for the natural convection of 
the ullage gas. The Nusselt number based on the tank radius as a function of Rayleigh number is plotted in Figure 
28. In general, one can see that for all the cases, heat transfer increases with Rayleigh number at a slow rate from 
low to intermediate Rayleigh numbers (<10 7 ). This is consistent with the variation of temperature stratification seen 
from the plots in Figures 26-27. With the increase of Rayleigh number, the temperature gradient increases near the 
gas-liquid interface, but the convection effect is rather negligible. For Ra > 1 0 7 , Benard-type cells start to develop 
near the interface, which leads to enhanced heat transfer. Just as for the case of liquid convection, the data are 
scattered, and there is no single correlation that can fit all the data. 

In an attempt to find a characteristic length that can correlate all heat transfer rates at different fill levels into a 
simple curve, we choose the vertical height of the ullage gas volume. It measures from the gas -liquid interface to the 
top of the dome. It is consistent with the characteristic length for the liquid convection. The results are shown in 
Figure 29. Here indeed, data collapse into a single curve. Using this length scale, the Nusselt number is almost 
constant for the value of Ra less than 10 7 . This is very much heat conduction regime. For Ra from 10 7 to 10 12 , 
Benard-type cells start to develop, and the curve follows a power of %, just as in the case of laminar heat transfer 
from a vertically heated plate using similarity solution. The general correlation can be expressed as: 


Nu l = 4.5 Ra L < 10 7 

Nu l = 0.08 (Ra L ) 1/4 10 7 < Ra L < 10 12 


(33) 


Figure 29 also shows the comparison of the presently developed correlation to the commonly used correlations 
built from 1960’s for vertically or horizontally heated flat plates. It is apparent that using the previous correlations 
can over-predict heat transfer rate by as much as one order of magnitude. 



Figure 28: Predicted Nusselt number as a function of Rayleigh number with tank radius as the length scale 
and comparison to the correlation for either vertically heated plate or horizontally heated plate. 
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Figure 29: The present correlation for natural convection of ullage gas inside a propellant tank using the new 
characteristic length and comparison to that of heated flat plate. 

VIII. Conclusion 

A fundamental first principles based approach has been undertaken to improve understanding of heat transfer in 
cryogenic tanks. The CFD code, CFD-ACE+, was first validated against available experimental data and commonly 
used correlations for natural convection along a vertically heated wall. Good agreements between the present 
predictions and experimental data were found for flows in laminar as well as turbulent regimes. 

With the validated CFD tool, the convective heat transfer between the tank wall and cryogenic propellant, and 
that between the tank wall and ullage gas have been simulated. In order to extract time independent heat transfer 
coefficient, the domain is divided into two parts, one for ullage, and one for liquid fluid. 

The present study shows that commonly used heat transfer correlations for either vertical or horizontal plate will 
over-predict heat transfer rate by as much as one order of magnitude. 

A characteristic length scale, defined as the wetted height L, has been defined that can correlate all heat transfer 
coefficients for different fill levels into a single curve. This curve can be used for the reduced heat transfer model 
analysis. 

In summary, the correlated heat transfer coefficient for the liquid propellant phase convection is: 


Nu l = 0.642(RaJ 1/6 

Nu L = 0 A 67 (Ra L ) 1 ' 4 

Nu l =0.00053(Ra L ) l/2 

For the ullage gas phase convection, it is: 

Nu l = 4.5 

Nu l =0.08 (Ra L ) U4 


Ra L <10' 


10 


10 7 < Ra L < 10 
10 10 < Ra L < 5x10 


(34) 


13 


Ra h < 10 7 

10 7 <Ra L <10 


12 


( 35 ) 
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